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Abstract 

We focus on the non-equilibrium two-bath spin-boson model, a toy model for examining quantum 
thermal transport in many-body open systems. Describing the dynamics within the NIBA equa- 
tions, applicable, e.g., in the strong system-bath coupling limit and/or at high temperatures, we 
derive expressions for the cumulant generating function in both the markovian and non-markovian 
limits by energy-resolving the quantum master equation of the subsystem. For a markovian bath, 
we readily demonstrate the validity of a steady-state heat exchange fluctuation theorem. In the 
non-markovian limit a "weaker" symmetry relation generally holds, a general outcome of mi- 
croreversibility. We discuss the reduction of this symmetry relation to the universal steady-state 
fluctuation theorem. Using the cumulant generating function, an analytic expression for the heat 
current is obtained. Our results establish the validity of the steady-state heat exchange fluctuation 
theorem in quantum systems with strong system-bath interactions. From the practical point of 
view, this study provides tools for exploring transport characteristics of the two-bath spin-boson 
model, a prototype for a nonlinear thermal conductor. 
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I. INTRODUCTION 



Impurity models were proved to be extremely useful in predicting many physical phe- 
nomena. The famous spin-boson model [1, 2], describing the dynamics a single charge on 
two states coupled to a dissipative bath, e.g., a solvent, exhibits rich phenomenology, in- 
cluding various phase transitions. Its applications range from charge transfer reactions in 
biological systems [3], photosynthesis [4], and the Kondo problem for magnetic impurities 
[5]. A variant of the model is the spin-fermion model, where a qubit (spin) interacts with 
one or more metallic environments [6-8]. These celebrated impurity models are appealing 
from various reasons. First, they enclose rich dynamical phenomenology, e.g., the Marcus 
theory [1] and the Kondo physics [5]. More recently, addressing molecular electronic ex- 
periments, such generic models were proved to be useful in predicting various aspects of 
molecular transport characteristics [9]. Secondly, they serve as a benchmark for develop- 
ing simulation techniques and approximation schemes, for describing the dynamics of open 
many-body systems [10, 11]. 

The traditional spin-boson (SB) model, considering an impurity-spin coupled to a single 
thermal reservoir, serves as a prototype model for exploring quantum dissipation problems 
[1]. The non-equilibrium version of this model, referring to the case where the spin (subsys- 
tem) is coupled to two thermal reservoirs, has been suggested as a toy model for exploring 
quantum transport phenomenology through an anharmonic nanojunction [12, 13]. In this 
case, the generic situation is one of a non-equilibrium steady-state, regardless of the initial 
preparation. We refer to this model as the "non-equilibrium spin-boson model" (NESB). 
Given the complex dissipative spin dynamics observed in the single-bath SB model [1], one 
expect its non-equilibrium extension to reveal tangled transport properties. Fundamental 
topics of interest are the scaling of the energy current with the spin-bath coupling strength, 
the role the reservoirs spectral function and the tunneling splitting on the subsystem dy- 
namics and the transport coefficients, and the onset of nonlinear current-temperature bias 
characteristics at strong interactions. 

The transport behavior of the unbiased (zero magnetic field) NESB model has been 
studied perturbatively, under the assumption of weak system-bath interactions, using master 
equation methods [13, 14]. While this scheme, providing simple analytic expressions, can 
capture some of the aspects of the energy transport process, the inherent weak system- 
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bath coupling assumption results in a resonance-sequential transport process where the two 
reservoirs separately excite and relax the subsystem. Exact numerical results can be obtained 
by following the Keldysh approach [15] or by using the complex machinery of the multilayer 
multiconfiguration Hartree approach [16]. 

In this paper we present an analytical study of the NESB model in the strong coupling 
limit and/or at high temperatures. In this limit a concerted action of the two baths takes 
place, where at each relaxation or excitation process both reservoirs contribute in a non- 
additive manner. This renders the master equation description complex, since the amount 
of energy transferred between the two baths is no longer in a one-to-one relationship with 
the number of spin flip events. The objective of our analysis is the cumulant generating 
function (CGF). With this at hand, one can derive analytic expressions for the transport 
coefficients: the current and its cumulants, exposing their dependence on the microscopic 
parameters. Furthermore, given the CGF, the validity of the steady-state heat exchange 
fluctuation theorem [17] can be established for anharmonic quantum models in the strong 
coupling limit. 

The fluctuation theorem (FT) for entropy production quantifies the probability of neg- 
ative entropy generation, measuring "second law violation" [18, 19]. Both transient and 
steady-state fluctuation theorems (SSFT) have been derived, where the former looks at non 
steady-state processes over a finite time t, and the latter measures entropy production in 
non-equilibrium steady-state systems over a long interval. In the context of heat exchange 
between two equilibrium reservoirs, v = L,R, the SSFT can be roughly stated as [17, 20] 

\n[Pt(+uj)/Vt(-uj)} = A0w. (1) 

Here Vt(w) denotes the probability distribution of the net heat transfer uj, from L to R, 
over the (long) interval t, with A/3 = Tr 1 — T L 1 as the difference between the inverse 
temperatures of the reservoirs. Extending the work and heat FT to the quantum domain 
has recently attracted significant attention [21, 22]. Specifically, a quantum exchange FT, 
for the transfer of energy between two reservoirs maintained at different temperatures, has 
been derived in Refs. [17, 23, 24] using projective measurements, and in Refs. [25, 26], based 
on the unraveling of the quantum master equation (QME). These derivations assume that 
the interaction between the two thermal baths is weak, and can be neglected with respect 
to overall energy changes. Using the Keldysh approach, an exact analysis was carried out 
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in Ref. [27]. However, it is valid only for harmonic systems. It is thus an open question 
whether a heat exchange FT is obeyed by an anharmonic quantum system strongly coupled 
to multiple reservoirs. 

Another subtle point is the role of non-markovian effects on the heat exchange SSFT 
and the current cumulants. In charge transfer problems, this topic has recently attracted 
significant interest [28-30]. The analogous problem, the reflection of non-markovian effects 
within the CGF in energy exchange scenarios has been considered for equilibrium systems 
in Ref. [31]. The Markov approximation is justified once the relaxation of the bath is fast, 
while the dynamics of the subsystem is slow. In this case the amount of energy transferred 
between the subsystem and the bath is pinned down with an arbitrary precision, as a strict 
energy conservation condition is enforced. However, once the assumption of markovianity 
is relaxed, when the dynamics of the baths degrees of freedom is on a comparable timescale 
with the subsystem evolution, energy-non-conserving processes on short time scales due to 
the energy-time uncertainty (when looking only at a subsystem) cannot be excluded. On 
this bath-decorrelation time scale, it is not obvious that the basic symmetry [Eq. (1)] still 
holds. 

Considering the NESB model in the strong interaction limit, allowing for non-markovian 
effects, it is our objective here to investigate its heat exchange properties: (i) To obtain the 
CGF and gain an explicit expression for the heat current, useful for understanding heat cur- 
rent characteristics for anharmonic-strongly coupled systems, (ii) Given the CGF, to derive 
the heat exchange SSFT. (iii) To understand the role of non-markovian (memory) effects 
on the onset of the SSFT. Our analysis makes use of the noninteracting-blip approximation 
(NIBA) [1] . This scheme can faithfully simulate the SB dynamics at strong system-bath in- 
teractions and/or at high temperatures in the Ohmic case. It is also exact for the unbiased 
case at weak damping. Under this approximation, the subsystem's dynamics is described 
within a time convolution quantum master equation. We unravel this dynamical equations 
into trajectories with a particular amount of net energy dissipated at each contact. In the 
markovian limit a heat exchange SSFT is verified. We also obtain the CGF, independent 
of the particular physical realization. In the non-markovian case a symmetry relation is 
recovered [22], reaching the universal SSFT once the observation time t [Eq. (1)] is much 
greater than the bath decorrelation time. 

The paper is organized as follows. In Sec. II, we describe our model and recall known 
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results for the spin-boson model in the strong coupling limit. Sec. Ill presents results for 
the CGF in the markovian limit, introducing the concepts and definitions that will become 
useful once the more involved non-markovian case is considered in Sec. IV. In Sec. V we 
conclude. 

II. MODEL AND DYNAMICS 

The non-equilibrium spin-boson Hamiltonian, comprising a spin subsystem coupled to 
two [y = L, R) independent phonon baths, maintained at a temperature T v , is described by 
the Hamiltonian (h = 1) 

H = Y az + f ° x + ° z S A ^>^> + h») + E w A f , A- (2) 

Here a x and a z are the Pauli matrices, uo is the energy gap between the spin levels, and 
A is the tunneling energy. Explicitly, in the two-state basis, a z = |1)(1| — |0)(0| and 
a x = |1)(0| + |0)(1|. Each reservoir includes a collection of uncoupled harmonic oscillators, 
bj (bj jV ) is the bosonic creation (annihilation) operator of the mode j in the v reservoir. 
The parameter accounts for the system-bath interaction strength. 

The transport characteristic of the non-equilibrium spin-boson model can be obtained 
exactly using numerical simulations [16]. Here, with the motivation to gain insight into 
the heat current characteristics, the behavior of the current cumulants, and the fluctuation 
symmetries we resort to approximations, allowing for analytical results. In particular, we 
employ the NIBA equations, valid at strong system-bath interactions or for high tempera- 
tures, assuming an Ohmic spectral density [1, 2]. The NIBA equations can be also obtained 
by applying the Born approximation with respect to the dressed tunneling elements [32, 33]. 
While this method has been originally derived for a spin coupled to a single bosonic reservoir, 
one can trivially generalize it to describe a multi-bath case. 

We begin by transforming the SB Hamiltonian (2) to the displaced bath-oscillators basis 
using the small polaron transformation [34], H p = WHU, U = e tCTzn / 2 , 

H P = y o z + | {a + e iQ + a^) + £ ^4 A- (3) 

where a± = ^(<j x ±ia y ), or cr + = |0)(1|, <r_ = |1)(0|, are the auxiliary Pauli matrices, 
n = J2u^i and ^ = 2 ?Ej (&],„ - b j,v)- Under the NIBA approximation [2, 32, 33], 
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generalized to the two-baths case, it can be shown that the spin polarization (cr z (t)) obeys 
a convolution-type master equation 

ft r t 



d{a z 



dt K s (t - r)(a z (r))dr - J K a (t - r)dr, (4) 

where the symmetric and antisymmetric kernels are given by 

K 8 (t) = A 2 e- Q ' {t) cos[Q"(f)] cos(u; t) 

K a (t) = A 2 e-«'W sin[g"(0] sm(u t). (5) 

The complex function Q(t) = *Yl u Qv{t)-, made of a real and imaginary components, Q u (t) = 
Q' u (t)+iQ"(t) : is defined by the correlation function e~ Q ^ = (e* n We~* n ^), with the thermal 
average performed over both reservoirs degrees of freedom. It is given by 

Q'lit) = [°° ^-sm(ut)duj, 
Jo ^ 

Q'u(t) = r ^^[l-cos(iot)][l + 2n u (io)]dio. (6) 

JO TtOJ 

Here J v {u) is the z/-bath spectral function, incorporating system-bath interactions 

J V ( U ) = 4*Y, X )A U - U i)- ( ? ) 

3 

In what follows we focus on the two-state population dynamics, therefore we rewrite Eq. (4) 
in terms of the states population 



dt 2 
+ 



f e -Q'(t-s) cos [ WQ ( t _ s ) _ Q»( t _ s )] Pl ( s )ds 
Jo 

^ / e- Q '^ cos[cu (t - s) + Q"(t - s)] Po (s)ds, 
^ Jo 



1 = p (t)+ Pl (t), (8) 

where (cr z (t)) = Pi(t) — P o(i). We explore next the heat transport characteristics under the 
NIBA approximation (i) assuming a markovian dynamics, and (ii) more generally, retracting 
to the non-markovian case, allowing for memory effects in the thermal baths. The non- 
markovian analysis can be reduced to the markovian description in the appropriate limit. 
For clarity, we have decided to first present here the (simple) markovian limit, then generalize 
the analysis and portray the non-markovian regime. This allows us to introduce the main 
concepts involved in the CGF derivation within a relatively simple setup. 
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III. MARKOVIAN LIMIT 



A. Population Dynamics 

A general analysis of counting statistics of a multi-state system connecting two non- 
equilibrium markovian reservoirs has been carried out in Ref. [35], based upon the NIB A 
equations. We use this scheme and derive here the CGF for the NESB model. In the 
markovian limit one assumes that the spin system slowly evolves in comparison to the reser- 
voirs evolution. Thus, we make the following two simplifications in the integro-differential 
equation (8): First, we replace the population, p n (s) by p n (t) (n = 0,1), supposing that 
the timescale over which the memory, represented by the integral, is important, is suffi- 
ciently short. Second, we extend the integral upper limit to infinity, assuming the integrand 
quickly dies out. Under these approximations, Eq. (8) reduces to a kinetic equation for the 
population dynamics, 

Pi = -kdPi(t) + k u p (t). (9) 
The rate constants are given as Fourier transforms of bath correlation functions, 

k d = C(co ), k u = C(-co ), (10) 

with 

/oo 
e^ ot C L (t)C R (t)dt. (11) 
-oo 

The ingredients of this correlation function are given in terms of the function Q v (t), defined 
in Eq. (6), 

C*(t) = | e -^W. (12) 

Using the convolution theorem, the transition rates C(±w ) can be rewritten as a convolution 
of the L-bath and i?-bath induced processes, 

1 f°° 

C(cj ) = — C L (cu -u)C R (u)du, (13) 

^ J-oo 

introducing the Fourier transform 

/oo 
e™C(t)dt. (14) 
-oo 
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These bath-specific microscopic rates satisfy the detailed balance relation, 

e wA \ (15) 



C v (-u) 

However, such a relation does not hold for the combined rate C(oo), ruling the dynamics. 
The QME (9) encloses complex physical processes as Eq. (13) draws nontrivial transfer 
rates. When the system decays it disposes the energy ooq into both reservoirs, cooperatively; 
the energy oo is dissipated into the R bath while the L bath gains (or contributes) the rest, 
oo — oo. Similarly, excitation of the system occurs through an L-R compound process. Since 
energy is dissipated or absorbed in such complex processes, energy "counting" is a nontrivial 
task, as reflected in the resolved master equation (16) discussed below. 

B. Cumulant Generating Function 

We construct next the cumulant generating function for the NESB model in the NIBA- 
markovian limit presented above. Following Ref. [26], we begin by defining the function 
Vt(n,u) as the probability that within the time t a total of energy oo has been transferred 
from the left bath to the right bath, while the spin is populating the n (n — 0, 1) state at 
time t. The time evolution of this quantity follows 

= -Vt(0,u) / —C R {-oo)C L {oo-oo )doo 

dt J_ OQ Ztc 

f°° 1 

+ / —C R (u-u)C L (u -(u-u))V t (l,u)du 

J-oo 27T 

= £° ^C R {u)C L {u -u)du 

f°° 1 

+ / —C R (u-u)C L ((u-u)-u )V t (0,u)du>, (16) 

for details see Appendix A. One can rationalize this equation as follows. Focusing for 
example on the dynamics of P*(l, oo), the first term in this rate equation describes the decay 
of this probability as the spin flips to the ground state and extra energy is dissipated to the 
R reservoir. The second term collects processes with an energy Co transferred to the R bath 
by the time t, with the spin occupying the ground state. At time t a spin flip takes place 
accompanied by an extra energy oo — Co dissipated to the R bath, completing the transfer of 
a total amount of energy oo to the right bath at t. 
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We note that in the present model we cannot adopt standard approaches for unraveling 
the reduced density matrix, directly dressing the interaction term in the Hamiltonian by the 
counting process [21]. This in because the two reservoirs affect the energy transfer process 
in a nonlinear way, thus counting system-bath interaction processes (as in a perturbation 
theory series) does not reveal here the actual amount of energy exchanged between the two 
reservoirs. 

We Fourier transform the above system of equations to obtain the characteristic function 
Z(x, t) for the energy counting field x, 

z (x,t))=[ J -~ . ] (17) 



It satisfies the differential equation 

d\Z( X ,t)) 



= -W(x)\Z(x,t)), (18) 



dt 

where the matrix W contains the following elements 

= l c( } _c*< A 
\-c(x) cm ) 

The diagonal terms were defined above, see Eq. (11). The nondiagonal terms are given by 
the integrals 

1 r°° 

C d/u ( X ) = — / C d/u (co)e luJX doo (20) 

with the components 

C d (cu) = C R (cu)C L (cu -u) 

C u (u) = C R (u)C L (-cu-cu ). (21) 

The cumulant generating function is formally defined as 

1 f°° 

G( X ) = lim - In / V t (co)e tu)X duj, (22) 
t 

where we introduced the short notation, V t {oo) = P t (0,u) + P t (l,u), the probability to 
transfer by the time t an energy u> from left to right, irrespective of the spin state. In the 
present case the CGF is expressed in terms of \Z) as 

G(x) = lim -\n(I\Z( X ,t)), (23) 

t— >oo t 



with (7 1 = (11 1, denoting a left vector of unity. It is practically given by the negative of the 
smallest eigenvalue of the matrix W. We diagonalize W and explicitly obtain the CGF in 
terms of the microscopic rates, 

G{x) = cm + c(-u ) + [{CM-c{-uj )f + AC d {x)c u {x)] l/2 (24) 

The heat current and its noise power can be readily derived, by taking the first and the 
second derivatives, respectively, of the CGF 



(J) 



(u) t _ dG{ X ) 



d(ix) 



{S) = i " W *=o (25) 

Here {u) t denotes the total energy u transferred from L to R by the (infinitely long) time 
t. Using the formal structure (24) one can show that the steady-state heat current, defined 
as positive when flowing left to right, obeys 

1 [°° 

( J) = — / udu [C r (uj)C l (u) - u)pi - Cr(-uj)C l (-ujo + u)p ] . (26) 

27r y_ 00 

This expression incorporates the steady-state populations 

Pl = C(-u )/(C(u ) + C(-u )), Po = C(u )/(C(u ) + C(-oj )). (27) 

For details see Appendix B. The result for the heat current agrees with the expression used 
ad-hoc in Refs. [12, 13]. It can be rationalized by viewing ujCr{u)Cl(ujq — ui)p\ as a spin 
relaxation process with the energy oj directed to the R bath and the amount of u — ui 
disposed into the L bath. Similarly, the second term describes energy loss from the R bath, 
where, combined with an energy influx from the L bath, results in the excitation of the 
spin system. It is significant to note that this expression has been achieved under relatively 
general conditions, for systems satisfying a markovian-NIBA approximation. The details 
of the Kernel K s / a {t) (e.g., the bath statistics) are not utilized in this derivation. Thus, 
it is valid for other systems following the structure (9)-(ll), e.g., the spin-fermion model 
[6, 7, 35]. Appendix B further details the derivation of the the second cumulant, the noise 
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power of the NESB junction, 

/DC pOO 
— lo 2 C r {uj)C l {ujq - u)duj + p —u 2 C r (-cj)C l (cj - Lu )du + 
■oo ^ J-oo 21T 

j | poo poo 

- 2 7v — \ i I uC R (-u)C L (u - uj Q )duj / ujC r (u)C l (uj - u)duj + 

C. Fluctuation Theorem 

We continue and confirm the validity of the SSFT in the NESB model, under the Markov 
approximation. This relation can be established by examining the symmetry of the CGF, Eq. 
(24) [21]. It is clear that it is sufficient to focus on the product term, V(x) = C d (x)C u (x), for 
resolving the symmetry of G(x)- Using the definitions (20)-(21) we therefore write (ignoring 
the 2ir prefactors) 

/oo poo 
e^C R (ou)C L {uj - u)du x / e^C R (cu)C L (-cu - u>)du>. (29) 
oo J —oo 

Shifting the argument x ~^ — x), A/3 — (3 R — /3l, it translates to 

/oo 
e-^e~^C R (u)C L (cu -u)duj 
-oo 

/oo 
e-™ x e-" Ap C R {u)C L {-u - u)doo. (30) 
-oo 

We now change variables, u — > —u, then use the detailed balance relation for C u (u), see 
Eq. (15). This transforms the first element in the RHS of the above equation to 

/oo 
e iux e^C R {-uj)C L {uj Q + u)du 
■oo 

/oo 
■oo 

= e huJ0 C u ( X )- (31) 
Similarly, the second element in the RHS of Eq. (30) reduces to 

/oo 
e ^^C R (-u)C L (u-u )du 
■oo 

/oo 
e^ x e^C R {u)e-^C L {uj Q - ^e-^-^dco 
■oo 

= e-^C d (x). (32) 
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Joining these two pieces we conclude that 

C d (x)C u ( X ) = C d (tA(3 - X )C u (tAf3 - x ). (33) 
Therefore, the CGF overall satisfies 

G(x) = G(iAP-x). (34) 

We are now in position to demonstrate the validity of a fluctuation relation for this non- 
equilibrium strongly coupled system. The probability to transfer the energy oj by the time 
t, from L to R is given by the inverse Fourier transform of Eq. (22), 

V t {u) = — / e tG(x) e-^ x dx- (35) 



2tt . 

Similarly, the quantity Vt(-ou) represents the probability that overall an energy u has been 
transmitted in the opposite direction, right to left, up to time t. Based on the symmetry of 
the CGF, Eq. (33), one readily concludes that 

l im li n Vt ^\ =cu Af3. (36) 

This expression describes a fluctuation relation for the non-equilibrium SB model, valid be- 
yond the weak-coupling approximation [25, 26]. Comparing this result to the weak coupling 
limit, described in Appendix C, we observe that formally these two expressions are iden- 
tical. However, one should note that in the strong coupling limit the energy variable oj is 
continuous, since multi-phonon processes in which part of the energy goes to L bath and 
part goes to the R baths, are allowed. In contrast, in the weak coupling limit energy transfer 
processes take place in integer units of the spin spacing, since this energy travels to either 
reservoirs separately. 



D. Examples and the Gaussian-Marcus limit 

We exemplify our results, and work out an expression for the heat current and the noise 
power for a specific case, the so called "Marcus" limit [3], assuming high temperatures 
T u > ujq and strong coupling. This limit is reached by performing a short time expansion of 
Q(t), [Eq. (6)] resulting in 

Q' v {t) = KT„t\ Q" v {t)=E»t. (37) 
12 



The reorganization energy E v r = £\ ^ 2 j,u/°°j incorporates system-bath interactions. It can 
be equivalently expressed in terms of the spectral density (7), E" — J ^^-dui. Using these 
expressions, the Fourier transform of the time dependent rates (11) and (12) can be resolved, 



7T 



cxp 



r - 1 - v 



CM ~ T\Ie?t r Ie?t l 



cxp 



(w - E^ - E, 



R\2 



4(T R E? + T L E$) 

Following Eq. (26), the average heat current can be analytically obtained 



(38) 



(J) 



2-itE^E^AT 



(2E r L T L + 2E?T R )i 



exp 



(E? + E?-u ) 2 
'4(Ef<T L + E?T R ) 



x f A 



(39) 



where Ja 



cxp 



;n7"L7'~iv: 1 1 . This result agrees with [12, 13]. The second 
moment of the current can be similarly calculated, however the expression is too cumbersome 
to be included. We present the behavior of the current and its noise power in Fig. 1. We 
observe nonlinear effects in the energy current, including the effect of negative differential 
conductance [12, 13]. The noise drops with increasing bias temperature. 

One can in principle seek to derive an analytic form for the probability distribution 
function Vt(uj), since the analytic structure of W is known: The diagonal terms are the 
rates, (38), the nondiagonal part is given by 



' {X) ~ 4 V E L T L + E R T R 



e 4(E L T L +E R T R ) ^40) 



Since the result is very complex, we retreat here to numerical simulations. We plug these 
expressions into the formal solution (24) and perform a numerical inverse Fourier transform 
[Eq. (35)], to obtain the distribution Vt(w). Fig. 2 demonstrates its shape at different 
times. The averaged current and noise agree with the values provided in Fig. I. We have 
also confirmed that the different curves indeed satisfy the SSFT (inset). 

We now go beyond the Marcus limit, and demonstrate the behavior of the current with 
the reorganization energy E r , quantifying system-bath coupling strength. We consider the 
unbiased case Uq = 0, and take an Ohmic spectral function J v {u) = ^ L we" w / 1Jc , identical 
for the two baths. Since we are interested in the energy current behavior for both weak 
and strong coupling strengths, we numerically calculate the elements in Eq. (26) using the 
definitions (6) and (14). The results are displayed in Fig. 3 showing a turnover behavior, 
where the current decays with increasing E r , at large values. It can be easily proved that 
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the weak coupling scheme can only produce a linear dependence of the current on the 
coupling strength [12, 13], see also Eq. (C18). The decaying behavior observed here is thus 
a fingerprint of the strong coupling limit. Similar results were reported in [16], using exact 
numerical simulations. Practically, this turnover behavior indicates that for maximizing the 
rate of energy transport in nanodevices one should work at the intermediate system-bath 
coupling limit. 




1 2 3 1 2 3 

T L" T R V T R 




FIG. 1: Energy current and noise power of the spin-boson model in the Marcus limit, Eqs. (38)- 
(40). Other parameters are T L = 5, T R = T L - AT, oj = 0.5, A/2 = 1, (a)-(b) E v r = 1, (c)-(d) 
E v r = 50. 




FIG. 2: The probability distribution V T (oS) at various times, t = 20 (full), t = 50 (dashed), t = 100 
(dotted). Other parameters are Tl = 5, Tr = 3, = 1, ujq = 0.5, A/2 = 1. Data was generated 
using the Marcus rates. The inset demonstrates for the same data the validity of the fluctuation 
theorem, with the slope of A/3 = 0.133, t = 20 (o), t = 50 (dotted), t = 100 (□). 
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FIG. 3: Energy current in the unbiased spin-boson model, Tl = 5, Tr = 3, ujq = 0, uj c = 10, 
A/2 = 1, numerically simulating (6) and the resulting current (26). The inset zooms on the weak 
coupling limit, displaying a linear dependency of the heat current on the reorganization energy. 

IV. NON-MARKOVIAN DYNAMICS 
A. Cumulant Generating Function 



We generalize here the results of the markovian analysis and derive the CGF for the 
non-markovian model introduced in Sec. II. A systematic formalism for analyzing non- 
markovian effects in charge transfer systems has been detailed in Refs. [28-30] . Here we 
adapt this scheme to describe energy transfer processes. Further, while a single counting 
field has been introduced in [28] , in the present model one should introduce two such fields, 
independently counting energy transmission at each contact. We begin the analysis by 
rewriting the equation of motion for the spin population (8) as 



dpiit) 
dt 



A 2 
-— ft 



7 iuo(t-s) -Q L (t-s) -Q R (t-s) 



Pi(s)ds 



A 2 f* 
+ — ft / r 



e uv >e ^ v 'e 
iu (t-s) e -Q L (s-t) e -Q R (s-t) po ^ ds ^ 



Po(t) +Pi(t), 



(41) 



where we made use of the symmetry properties of the Q(t) function, Q'(t — s) = Q'(s — t) 
and Q"{t — s) = —Q"(s — t), see the explicit expressions in Eq. (6). Here ft denotes the real 
part. In the next step, we use the Fourier transform relation 



A 



At) 



1 
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e- luJt C u (cu)du, 



(42) 
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and write 



t poo poo 



= ^3? / ds Pl (s)e iulo{t - s) I dwiCiMe-*" 1 **-') / d^C l; U 2 )< ; ~ 



2vr 2 7 

1 pt poo poo 

+ ^ J ds Po (s)e-^ s ^ J du 1 C L (u 1 )e- iU1 ^ J du 2 C R (uj 2 )e 



oo J — oo 

t poo poo 

iui2(s—t) 



(43) 



We now energy-resolve this equation, p n (t) = J™ dui doj R Vt(n, u>l, ujr), looking for the 
probability Vt(n,ujL,ou R ) at the time t the spin occupies the n (n — 0, 1) state, an overall 
energy u R has been transferred to the R bath and oo L has been transferred to the left bath. 
Note that unlike the Markov case, we separately count the energy dissipated at each bath. 
This probability satisfies the differential equation 



dt 2tt 2 

^ ft poo poo 

-} 



i pt poo poo 

— / dsP 9 {l,u L ,u R ) / / d^d^CLMC^^e 1 ^-^-^'^] 

K JO J -oo J -oo 



2vr 



i pt poo poo 

— ds / c^dc^O^a^CflK^ 

K JO J -oo J -oo 

(44) 



where we used the fact that C„(u)) is a real function. An analogous equation can be written 
for the time evolution of the probability V t (0, u>l, u> r ). We now introduce two counting fields 
Xl an d XRj f° r each reservoir, and Fourier transform the above equation with respect to these 
two fields. Further, we Laplace transform the resulting equation, H(z) = J °° e~ zt h{t)dt. 
Utilizing Fourier transform and Laplace transform convolution relations, Eq. (44) reduces 
to a linear equation 

\Z(xl, XR, z)) = 7t77~ r ( z \ Z (Xl, Xr, *»*->oo ' ( 45 ) 

z - W(xl,Xr,z) 

with the initial value theorem invoked, h(t = 0) = lim z _ >00 zH(z). The vector \Z) is defined 
by 

, 7( » I Jo°° dte ~ Zt Too Vt(0, u L , u; R )e™e™ duj L du R \ 

\Z{Xl,XR,z)) = ,,00, x ■ ( 46 ) 

^ / °° dte~ zt f_°^ V t (l, u L , u; R )e™e l ™ du L du R ) 
and the kernel W represents the matrix 

T iw v / -l + (z) a-( X L,XR,z) 

W{ X l,Xr,z) = I x x | (47) 

\a + {XL,XR,z) -7 {z) 
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with the elements 



oo poo 



a + ( X L,XR,z) = 7^2 I I e luJlXL e l0J2XR C R (uj2)C L (uJi) ~ d^idc^a 

2vr 2 > /_ 00 > /_ 00 ^ 2 + (w + wi + w 2 ) 2 

= / / e^V™^^)^)^^ rjdwidwa 

1 f 00 f™ z 

7 + (X> = / Cr(w2)C l (wi)^— -^du^ 

2?H > /_ 00 > /_ 00 ^ 2 + (w + wi + w 2 ) 2 

(48) 

We are interested in the total probability, to occupy either states, 

V t (oj L ,u R ) = ^ P t (n,w L ,c^). (49) 

n=0,l 

We express it in terms of the characteristic function e s ^- XL ' XH,t \ 

/oo 
V t (u>L, oj R )e lulLXL e lulRXR duj L duj R . (50) 
-00 

This expression generalizes Eq. (22) to the non-markovian case. It can be formally expressed 
by an inverse Laplace transform of Eq. (45) [28], 



e S{xL ' XR ' t] = 7T- / dze«(I\ w r(*£(XL,X*, *)>),_►«, ■ (51) 

Here a is a real number, larger than the real part of all the singularities of the integrand. 
Equation (51) is a formal result. In practice, it is evaluated as follows: First, we note 
that the stationary solution (assumed to be unique) of Eq. (44) is given by the eigenvector 
corresponding to the zero eigenvalue of W, 

W(xl = 0,Xr = 0,z = 0)\Zss}=0. (52) 

Furthermore, as a result of the normalization and conservation of the total spin probabilities, 
an eigenvalue of W satisfies \q(xl = 0, xr — 0, z) = 0, for all z [28]. This can be directly 
verified in our case, Eq. (48): The element a + (a~) reduces to 7+ (7") when xl,r = 0, and 
a zero eigenvalue sustains, irrespective of the value of z. At finite value for the counting 
fields an eigenvalue A (xl? XRi z ) adiabatically develops from this zero eigenvalue, with small 
Xl, Xr an d z. The long time behavior of the characteristic function is therefore determined 
by the pole structure of (z — A (xi, X2, z))" 1 close to zero. This pole z (xl, Xr) solves 

z = X (Xl,Xr,z ), (53) 
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and it should reduce to z (xl = 0, xr — 0) = 0, describing the stationary state. Since 
all other singularities have larger negative real parts, this pole determines the long time 
behavior of the characteristic function as 

e s(xL,xn,t) _^ XRj Zo yo(xL,x R )t_ ( 54 ) 

In the markovian limit Ao does not depend on the z variable, thus trivially z (xl,Xr) — 

Concluding, the scheme to obtain the CGF proceeds as follows [28]: (i) We obtain 
^o{Xl, XRi z )i the eigenvalue of W(xl, Xr, z ) that adiabatically develops from the zero (sta- 
tionary) eigenvalue, (ii) We solve Eq. (53) and obtain the pole z (xl,Xr)- ( m ) We identify 
the CGF, the analog of Eq. (22), by the pole, 

G(xl, Xr) = z (xl, Xr)- (55) 
Back to (47), we resolve the eigenvalue 

\ i \ 7 + + 7~ . V(7 + -7") 2 + 4a+a- 

MXl,Xr,z) = g + 2 ' ^ ' 

satisfying \ (xl = ®,Xr = 0,2;) = for all z [28]. The elements of Eq. (56) all depend on 
the variable z, further depend on the counting fields. In principle, we should now solve 
Eq. (53) in order to gain the CGF, thus the current and its cumulants. 



B. Heat Current 

In the long time limit the combination of Eqs. (51) and (54) leads to 

z (xl,Xr)^ 7 In/ Vt(u L ,u R )e™™e™«x*du> L du> R . (57) 

It is argued in Ref. [28] that the current and its cumulants can be obtained directly through 
the analysis of A itself, Taylor expanded around z = 0, xl = 0, and xr — 0, 

^xM=J2 {iX f iiX f Z >^\ (58) 

n\ ml l\ 

n,m,l 

with 

c (n,m,l) = d^ L) d^ R) d l MXL, XR, z)\ Xl , Xr ,^ - (59) 
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The thermal current, calculated by counting energy flow at the L contact, is given by 



t 



.(1,0,0) 



Similarly, the current detected at the R contact satisfies 



(Jr) 



("tit 



.(0,1,0) 



(60) 



(61) 



or explicitly 



(Jr) = (7 + + 7-)- 1 



da + 

OL irr. r + OL 



da 



d(ix 



R) 



d(ix 



R) 



(62) 



X£,Xfl,z=0 

It can be easily verified that this two quantities are identical (with opposite sign), and 
equivalent to the markovian result (B4). It is thus significant to note that our formalism 
provides a general expression for the energy current, for many-body systems satisfying the 
dynamics (4), irrespective of the details of the thermal reservoirs. One can similarly calculate 
high order cumulants, by evaluating high order c terms [28]. 



C. Fluctuation Theorem 



The assumption of no memory enforces a strict energy conservation condition for processes 
between the system and the environments. In contrast, in the non-markovian regime there 
is no such an energy-conservation statement, thus it is not obvious that the general FT 
symmetry (34) still holds for any time interval t [22, 23]. We now prove that the eigenvalue 
^o(Xl, Xr-i z ) satisfies the symmetry relation 

Vxl, Xr, z ) = V^l - Xl, i/3r - Xr, z). (63) 

Only in the markovian limit the symmetry is given in terms of the affinity, as Xo(x) = 
\ (iA(3 — x)- Since the counting fields and z are independent variables, these symmetry 
relations translate into the analogous relations for the CGF itself, zq(xl,Xr)- This result 
exemplifies that while microreversibility is sufficient for deriving the basic symmetry relation 
(63), the SSFT holds only under more restrictive conditions, dictated here by the bath 
relaxation timescale [21, 23]. 

The symmetry of \o(xl, Xr, z )-> thus the symmetry of the CGF, is coded in the product 
of terms that depend on the counting fields, T>(xl, Xr-, z ) = ct + a~, see Eq. (56). We can 
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readily confirm that 

a + (il3 L - Xl, i/3r - Xr, z) = 

1 roo roc 

— / / e-^e-^e-^e-^CRMC^)^—. ■ ■ ^duu^ 

2vr 2 > /_ 00 J_ OQ z 2 + (u + uji + u 2 ) 2 

= ^(xl,Xr,z). (64) 

This result is obtained by changing variables, u± — > —Ui and u 2 — > —u 2 , then utilizing the 
detailed balance relation, C v {— cu) = C v (w)e~^ . Similarly, it can be proved that 

a~(i/9 L - Xl, ipR - Xr, z) = a + (x L , Xr, z ). (65) 

As a result, the symmetry relation (63) is confirmed, and the CGF, reached by solving Eq. 
(56), similarly satisfies 

z (xl, Xr) = zo{i(3 L - Xl, i/3r - Xr)- (66) 

The probability itself is given by an inverse Fourier transform [Eq. (50)] with respect to 
both fields, 

-1 poo poo 

V t (u L ,cu R ) = —— / / e^^e-^e'^^dxLdxR- (67) 

W J-ooJ~oo 

Based on the symmetry of the CGF, it can be readily proved that in the long time limit the 
following "basic" fluctuation relation holds [22] 

V t (u L ,u R ) = e p LULe p RWR , Q8) 
V t {-oo L ,-u R ) 

The "standard" fluctuation theorem, expressed in terms of the affinity A/3 = (3 R — (3l is 
regained when the kernel W reduces to the markovian result; The two counting fields then 
trivially count the same amount of energy, oj l = —uj r . This can be explicitly shown by 
evaluating the elements a + and a~ in the markovian limit z — 0. We find that 



a 



1 r°° 

(Xl, Xr,z = 0) = — / e-^ + ^e^C R (u 2 )C L (-u - oo 2 )du 2 
1 f°° 

(Xl, XR,z = 0) = — / e^-^e^C R (u 2 )C L (uo - u 2 )du 2 
2tt J_ 00 



(69) 



We now define a new counting field, X = Xr ~ Xl, and immediately verify that the product 
of these two objects, V, satisfies 

V(x)=V(i/\(3-x)- (70) 
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This directly implies on the same symmetry for the markovian CGF, 

Mxr-Xl) = X (iA/3 -xr + Xl), (71) 

leading to the standard heat exchange SSFT, Eq. (36). 

One should note that Eq. (54) has already involved the assumption of long times, such 
that only one eigenvalue of W, with the smallest (absolute) real value, dictates the dynam- 
ics. The z dependence in Eq. (48) thus manifests itself when the bath decorrelation time 
is long, comparable with the inverse relaxation rates of the subsystem. This observation 
establishes the regime of validity of the SSFT, Eq. (34). It holds when the interval t is long, 
beyond the bath memory time. We recall that for strictly harmonic systems one directly 
obtains the SSFT [27], without any reference to the bath characteristic timescale. This 
is because in coherent systems the reservoirs only serve as a source for excitations, which 
then elastically cross the impurity. In contrast, in the present model inelastic bath-induced 
processes are involved in the energy transfer process, making the bath decorrelation time a 
relevant parameter for the dynamics. 



V. CONCLUSIONS 



We presented here a scheme for obtaining the CGF, thus the current and its moments 
for the non-equilibrium spin-boson model, an eminent many-body impurity model. A heat 
exchange SSFT was established for quantum systems incorporating strong system-bath in- 
teractions and anharmonic effects. Our derivation relays on the NIBA equations, originally 
developed for the equilibrium spin-boson model, generalized to describe the dynamics of a 
spin impurity coupled to multiple thermal reservoirs. Our study provides closed expressions 
for the CGF, useful for deriving the distribution of heat fluctuations, the averaged current 
and the thermal noise power. We also showed explicitly that the timescale controlling the 
onset of the SSFT is the decorrelation time of the reservoirs. Future work will be devoted 
to generalizing our study to systems showing coherence effects, either using path integral 
formulation, or quantum master equation methods. Exploring the analogous dynamics for 
a fermionic system under voltage and temperature biases will be the topic of future studies. 
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Appendix A: Derivation of EOM for the resolved probability 

The equation of motion for the resolved probabilities in the markovian limit, Eq. (16), are 
explained here, based on the population dynamics (9). For clarity, we include this equation 
again, 

= -pi(t)C(u )+ Po (t)C(-ujo). (Al) 

The population of each state at time t can be expressed in terms of the resolved probability 
Vt(n,u), n = 0, 1, that within the time t a total of energy u has been transferred from the 
left bath to the right bath, while the spin is populating the n (n — 0, 1) state at time t, 

/oo 
V t (l,u)du>, 
■oo 

/oo 
Pt(0,u)du. (A2) 
■oo 

Plugging these integrals in the dynamical equation (Al), it becomes (ignoring 1/2% factors 
for simplicity) 

d Z* 00 /' oc /»oo 

— / Vt(l,uj)du = - C L (u - uijC^uijdu-v x / V t (l, uj 2 )duj2 

tit J — oo J — oo J — oo 

/oo ;>oo 
C L (-u - uJCnMdun x Vt(0,u 2 )du 2 . (A3) 
■oo J — oo 

We now equate identical energy terms, thus get the resolved dynamics (16) 

^ = -Vt{l,U,) J" CLfa-LjjCRMdU! 

/oo 
V t (0, wi)Cl(-w -u + ui)C r (w - wi)dwi. (A4) 
-oo 

For further validating this equation, we attempt to recover Eq. (Al) by integrating this 
equation over frequency. The first term in Eq. (A4) trivially reduces to the first term in 
(A3). The second term in Eq. (A3) is restored following a variable change, 

/oo roc 
du / V t (0,u>i)C L (-L} - u + u 1 )C R (u -u 1 )du 1 = 
oo J — oo 

/OO POO 
duj 2 C L (-ujo - u 2 )C R (u 2 ) / duV t (0,uj - u 2 ) = C(-u )p (t) (A5) 
-oo J — oo 
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Appendix B: Derivation of the current and its noise power in the markovian limit 



We derive here the steady-state heat current, Eq. (26), and the noise power Eq. (28). 
We begin by solving the kinetic equations (9) in the long-time limit. The steady-state 
populations are given by 



hi. 



Pi 

Po = 



C(-ouo) 



K + k d C(-u ) + C(u ) 
k d C(u ) 



(Bl) 



k u + k d C(-uo) + C(u ) 
We now study the first derivative of the CGF, Eq. (24), with respect to the counting field, 

dG( X ) 



(J) 



d(ix) 



d%x 



-c u ( x ) + 



d%x 



x [(C(u )-C(u )) 2 + 4C d (x)C u ( X )] 



-1/2 



(B2) 



Note that C d (x = 0) = C(u> ) and C u (x = 0) = C(-u ), a direct result of Eqs. (20) and 
(21). We identify the second term in the expression above by the sum (C(uj ) + C(— a; )) _1 . 
The partial derivatives are given by (ignoring (2ir) 1 factors for simplicity) 

dC u { X ) 



d(ix) 
dC d ( X ) 



d(i X ) 



/OO l>OQ 
luC 'r{uj)C — ooQjdu — — I ujC 'r{—uj)C 'l(u — u )du 
OO J — OO 

ojCr(uj)Cl(ujq — u)du. (B3) 



'OO 
OO 



Plugging these terms into Eq. (B2) we find the current 

1 



(J) = 



C(u ) + C(-uo) 

POO 

C(u) ) / ujC 'r(-u)C ' l (uj - ujojduj 



C(-u ) / <juC r (u)C l (uq - u)du 



(B4) 



It is significant to note that this expression stays intact for non- markovian systems [28]. 
Next we adopt the steady-state population (Bl) and simplify the result, 



/oo poo 
LdC R (cj)C L (u - uj)duj - po ujC r (-uj)C l {uj - cj )du, 
-OO J —oo 



(B5) 



which is precisely Eq. (26). We now verify that (J (AT)) = — (J(— AT)) for a spatially sym- 
metric system. Upon exchange of the temperature polarity, the above expression becomes 



/oo roc 
ujC l (uj)C r (uj q - uj)du -p uC l (-uj)C r (uj - u )du. 
-oo J —oo 



(B6) 
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We change variables, oo — ujq = —x, and get 

/oo 
(u - x)C l (cj - x)C R {x)dx 
-oo 

This expression can be organized as 
(J(-AT)) = u 



r 

-Po 

J — c 



(u - x)C L (x - u} )C R (-x)dx 

(B7) 



pi / C l (cj - x)C R {x)dx -p C L (x - u )C R (-x)dx 

J — oo J — oo 

/oo /»oo 
xC L (w - x)C R (x)dx -p / xC L (x - u )C R (-x)dx 
-oo J —oo 



• (B8) 



Since the first line fade away once combining the definition (11) and the steady-state pop- 
ulation (27), we establish the odd symmetry for the current with AT. The noise power is 
formally given by 



(S) = 



d 2 G( X ) 



d(i X Y 



d2Cd{x) C u ( X ) + d2CU{x) C d ( X ) + 2 d ° U{x) dCd{x) 
d{ix) 2 d{ix) 2 di X di X 

x [(C7(wb) - C(-coo)) 2 + AC d {x)C u {x)] ^ 



dCd(x) -c u ( x ) + dCU(x) 



c d {x) 



[(C(u )-C(-u )) 2 + 4C d (x)C u (x)} 



-3/2 



d(ix) ' ' d (ix) 
Using the explicit expressions for the correlations we reduce it to 



(B9) 

x=o 



(S) 



/OO /'OO 
u 2 C R (u)C L (u - u)duj + p uj 2 C r (-uj)C l (uj - uj )duj + 
-oo J —oo 

POO 

cuC R (— uj)Cl{uj — Lu )dcu / ooC R {uo)C l{uoq — u)dw + 



- 2 



1 

cm + c(-uj ) 

1 .-2 



CM + C7(-wo; 



(BIO) 



Appendix C: The spin-boson model in the weak coupling limit 



We study here the counting statistics of the unbiased (ojq = 0) spin-boson model, and 
verify the validity of the SSFT in this case, both under the Born-Markov Approximation 
[25, 26]. Our starting point is the SB Hamiltonian [Eq. (2)]. We take u = and apply a 
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unitary transformation 

U^a z U = <t x , U ] a x U = a z (CI) 

with U = -^(cx + CTz), to obtain the transformed Hamiltonian Hw = U^HU, 

Hw = Hq + Hi + Hb 
H = —(r z ; Hi = a x ^2 *j>( 6 J> + 6 j>) 

H B = H^Ys^A^- (C2) 

^ j 

Note that the subsystem energy gap is now given by A, rather than wo as in the original 
spin-boson description [Eq. (2)]. This form is a convenient starting point for a perturbation- 
theory calculation, assuming weak system-bath coupling. We outline next the principles of 
this standard approach [1]. Details, for the two-bath scenario, can be found in Ref. [13]. 
We begin with an equation of motion for the total density matrix p(t) in the interaction 
representation, 

p = -t[Hi(t),p(t)]. (C3) 

The operators are given by 0(t) = e ^ Ho+HB ' lt Oe^ H " +HB ' >t . We now make the following 
assumptions: (i) At the initial time the reservoirs are (separately) maintained in thermal 
equilibrium, isolated from the subsystem, and (ii) the spin and the baths are weakly coupled, 
allowing for a weak-coupling expansion with respect to the interaction term in Eq. (C2). 
This results in 

Ps(t) = -[ dsTr[Hi(t),[Hi( S ),ps(s)®p B }}, (C4) 
Jo 

which is a non-markovian equation of motion. Here ps(t) = Tr[p(t)] is the reduced density 
matrix of the system; the trace is performed over the two reservoirs. The bath density 
matrix is a product state, ps = Pl <8> Pr, of the two canonical density matrices p u = 
e~ Hv l T " j r \^ v \e~ H "l Tv \. In our model (C2) the population dynamics becomes decoupled from 
the coherences dynamics. It obeys 

p 1 = -2ul e lA{t - s) g(t -s) Pl (s)ds + 23? / e - lA(t - s) g{t - s)p {s)ds. (C5) 
Jo Jo 
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Here p n (t) = (ps(t)) n>n (n = 0, 1) and g(r) = J2 v 9u(t), with 

9u{r) = ^ [n u (u)e^ T + (n v (u) + l)e-^ T ] du. (C6) 

The bath spectral function is given by J u (uj) = 47r^- A| ^^(u; — Wj) and the function 
n u (u) = [e 13 ^ — denotes the Bose-Einstein distribution. We now make the marko- 
vian approximation, assuming that bath correlations decay on a time scale shorter than the 
subsystem characteristic timescale. This converts Eq. (C5) into a kinetic-Master equation, 



V 0->1 



Pi = -pi Yl +poJ2k 

V V 

Pi(*)+Po(*) = l, (C7) 

where the Fermi-golden rule transition rates are evaluated at the subsystem energy gap A, 
satisfying 

k^ = r„(AH(A), kUo = r,(A)[l + n,(A)]. (C8) 

The rate r„(u;) = 2rr ^ . A^5(wj — w) denotes the temperature independent part of the 
relaxation rate. The dynamics (C7) describes spin flip processes accompanied by an energy 
transfer at the amount of A to either the left or the right reservoirs. 

We proceed and derive the cumulant generating function in the present weak coupling 
limit following [26]. We begin by defining V t (n, qA) as the probability that within the time 
t a total energy gA has been transferred from the left bath to the right bath, while the 
spin is populating the n (n = 0, 1) state at time t. Note that q here is an integer, since 
energy is transferred here in discrete quanta of A, between the two baths. In other words, 
whenever the spin flips, the spin gap A is dissipated or absorbed at either the left or the 
right reservoir. Thus, 



dVti0 dt qA) = -Vt^qAXk^ + k^+Vtil^q-^k^ + VtihqA)^ 

= -V t (l,qA)(kL o + k^ o )+V t (0 : (q + l)A)k^ 1 + V t (0,qA)k^ v (C9) 



dt 

We Fourier transform these equations with the counting field x t° obtain the characteristic 
function, 

\J2 q V t (l,qA)e^ 
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satisfying a first order differential equation, 

d\Z( X ,t)) 



dt 



-W(x)\Z( X ,t)), 



(Cll) 



with the matrix 



W 



J- 

— h L — h R p"'X a 

^0-5>l ^0-5>l e 



(C12) 



The CGF is given by the negative of the smallest eigenvalues of this matrix, 



Gix) = 



-A + ^A 2 + AB(x) 



(C13) 



The coefficients are defined as 



A = r L [l + 2n L (A)]+r fl [l + 2n fl (A)], 
B( X ) = T L T R n L (A)n R (A) [(e~^ A - l)e feA + (e^ A - l)e feA ] . 



(C14) 



For brevity, we have discarded the direct dependence of the rates on frequency, T(A). It 
can be easily verified that the cumulant generating function satisfies the symmetry G(x) = 
G(iA/3 — x) w ith A/3 = /3 R — j3 L . This symmetry can be translated into the fluctuation 
relation at long time t [21, 26], 

V t (qA) 



= e 



(C15) 



V t (-qA) 

Comparing this result to the strong coupling expression (36), we note that in the unbiased- 
weak coupling limit the discrete energy qA replaces the continuous variable u, since the 
reservoirs here accept or contribute energy in quanta of the spin spacing A. Finally, the 
current and the noise power are given by 

1 dB 



(J) 



A dix 



The elements in this expression are 
dB 



d 2 B 



d(ix) 2 A 2 \dix 



dB 



(C16) 



d(ix) 

d 2 B 



x=o 



AT L T R [n L (A)-n R (A)}, 
: -A 2 T L T R [n L (-A)n R (A) + n R (-A)n L (A)]. 



(C17) 
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We find that in this weak coupling limit the current satisfies 



/ t\ a T L T R [n L (A)-n R (A)] 

{J) r L [l + 2n L (A)] + r JI [l + 2n Jl (A)]- 1 J 



This result agrees with previous studies [12]. 
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